args <- commandArgs(TRUE)
if(length(args) != 4){
	cat ("Usage:\n\tRscript base.R <Base_distributions_by_read_position_1.txt> <Base_distributions_by_read_position_2.txt> <output raw png file> <output clean png file>\n")
	q(save="no")
}
rt1_1 <- read.table(pipe(paste("cut -f2-6 ",args[[1]])),skip=1)
rt1_2 <- read.table(pipe(paste("cut -f7-11 ",args[[1]])),skip=1)
rt2_1 <- read.table(pipe(paste("cut -f2-6 ",args[[2]])),skip=1)
rt2_2 <- read.table(pipe(paste("cut -f7-11 ",args[[2]])),skip=1)
rawa <- c(as.numeric(sub("%","",rt1_1$V1)),as.numeric(sub("%","",rt2_1$V1)))
rawc <- c(as.numeric(sub("%","",rt1_1$V2)),as.numeric(sub("%","",rt2_1$V2)))
rawg <- c(as.numeric(sub("%","",rt1_1$V3)),as.numeric(sub("%","",rt2_1$V3)))
rawt <- c(as.numeric(sub("%","",rt1_1$V4)),as.numeric(sub("%","",rt2_1$V4)))
rawn <- c(as.numeric(sub("%","",rt1_1$V5)),as.numeric(sub("%","",rt2_1$V5)))
cleana <- c(as.numeric(sub("%","",rt1_2$V1)),as.numeric(sub("%","",rt2_2$V1)))
cleanc <- c(as.numeric(sub("%","",rt1_2$V2)),as.numeric(sub("%","",rt2_2$V2)))
cleang <- c(as.numeric(sub("%","",rt1_2$V3)),as.numeric(sub("%","",rt2_2$V3)))
cleant <- c(as.numeric(sub("%","",rt1_2$V4)),as.numeric(sub("%","",rt2_2$V4)))
cleann <- c(as.numeric(sub("%","",rt1_2$V5)),as.numeric(sub("%","",rt2_2$V5)))
pos1 <- seq(1,length(rawa))
pos2 <- seq(1,length(cleana))
png(file=args[[3]],height=360,width=576,antialias='default')
plot(pos1,rawa,main="Base percentage composition along reads",type="l",col="red",lwd=2, bty="o",xaxt="n",yaxt="n",xlab="",ylab="",xaxs="i",ylim=c(0,50))
lines(pos1,rawc,lty=2,col="green",lwd=2)
lines(pos1,rawg,lty=3,col="blue",lwd=2)
lines(pos1,rawt,lty=6,col="magenta",lwd=2)
lines(pos1,rawn,lty=4,col="cyan",lwd=2)
lines(c(length(rt1_1$V1),length(rt1_1$V1)),c(0,50),col="blue",lty=4,lwd=2)
axis(side=1, seq(0,max(pos1),20), tcl=0.2, labels=FALSE)
axis(side=2, seq(0,50,10), tcl=0.2, labels=FALSE)
axis(side=3, seq(0,max(pos1),20), tcl=0.2, labels=FALSE)
axis(side=4, seq(0,50,10), tcl=0.2, labels=FALSE)
mtext(seq(0,max(pos1),20),side=1,las=1,at=seq(0,max(pos1),20),line=0.3,cex=1.5)
mtext(seq(0,50,10),side=2,las=1,at=seq(0,50,10),line=0.3,cex=1.5)
mtext("Position along reads",side=1, line=2, at=max(pos1)/2, cex=1.5)
mtext("Percent",side=2, line=2.5, at=25, cex=1.5)
legend("topright",lwd=2,lty=c(1,2,3,6,4),col=c("red","green","blue","magenta","cyan"),bty="n",legend=c("A","C","G","T","N"),cex=1,inset=0.05)
d <- dev.off()
png(file=args[[4]],height=360,width=576,antialias='default')
plot(pos2,cleana,main="Base percentage composition along reads",type="l",col="red",lwd=2, bty="o",xaxt="n",yaxt="n",xlab="",ylab="",xaxs="i",ylim=c(0, 50))
lines(pos2,cleanc,lty=2,col="green",lwd=2)
lines(pos2,cleang,lty=3,col="blue",lwd=2)
lines(pos2,cleant,lty=6,col="magenta",lwd=2)
lines(pos2,cleann,lty=4,col="cyan",lwd=2)
lines(c(length(rt1_2$V1),length(rt1_2$V1)),c(0,50),col="blue",lty=4,lwd=2)
axis(side=1, seq(0,max(pos2),20), tcl=0.2, labels=FALSE)
axis(side=2, seq(0,50,10), tcl=0.2, labels=FALSE)
axis(side=3, seq(0,max(pos2),20), tcl=0.2, labels=FALSE)
axis(side=4, seq(0,50,10), tcl=0.2, labels=FALSE)
mtext(seq(0,max(pos2),20),side=1,las=1,at=seq(0,max(pos2),20),line=0.3,cex=1.5)
mtext(seq(0,50,10),side=2,las=1,at=seq(0,50,10),line=0.3,cex=1.5)
mtext("Position along reads",side=1, line=2, at=max(pos2)/2, cex=1.5)
mtext("Percent",side=2, line=2.5, at=25, cex=1.5)
legend("topright",lwd=2,lty=c(1,2,3,6,4),col=c("red","green","blue","magenta","cyan"),bty="n",legend=c("A","C","G","T","N"),cex=1,inset=0.05)
d <- dev.off()
